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Abstract 

The dynamics of generalized Lotka-Volterra systems is studied by theoretical techniques and 
computer simulations. These systems describe the time evolution of the wealth distribution of 
individuals in a society, as well as of the market values of firms in the stock market. The individ- 
ual wealths or market values are given by a set of time dependent variables W{, i = 1, ...N. The 
equations include a stochastic autocatalytic term (representing investments), a drift term (repre- 
senting social security payments) and a time dependent saturation term (due to the finite size of 
the economy). The w^s turn out to exhibit a power-law distribution of the form P(w) ~ w~ 1 ~ a . It 
is shown analytically that the exponent a can be expressed as a function of one parameter, which 
is the ratio between the constant drift component (social security) and the fluctuating component 
(investments). This result provides a link between the lower and upper cutoffs of this distribution, 
namely between the resources available to the poorest and those available to the richest in a given 
society. The value of a is found to be insensitive to variations in the saturation term, that represent 
the expansion or contraction of the economy. The results are of much relevance to empirical stud- 
ies that show that the distribution of the individual wealth in different countries during different 
periods in the 20th century has followed a power-law distribution with 1 < a < 2. 

PACS numbers: PACS: 05.40.Fb,05.70.Ln,02.50.-r 
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I. INTRODUCTION 



In recent years there has been considerable interest in the collection and analysis of large 
volumes of economic data. Such data includes the distributions of the income and wealth 
of individuals |1|, 0, [| f|, ||, |6|, 0, |^], the market values of publicly traded companies as well 
as their short and long term fluctuations 0, 10, 11, T2[. A common observation is that 



distributions of economic data exhibit a power-law behavior of the form 

P{w)^ W - 1 - a . (1) 

where the variable w represents the wealth of an individual or market value of a company 
and a is the exponent that provides the best fit to the empirical data. Empirical studies 
show that the distribution the wealth of individuals in different countries follows the power- 
law behavior described by Eq. ([[]), with 1 < a < 2 jl], @, 0, §|, |§, §| . These results stimulated 
theoretical studies in attempt to construct models that reproduce the power-law behavior 
and predict the value of a |]| 0, |TJ, [H| [T7|, |T? 



In this paper we study a stochastic dynamical model, based on the Lotka-Volterra system 
that gives rise to the power-law distribution of Eq. (|l|). The model consists of coupled 
dynamic equations which describe the discrete time evolution of the basic system components 
Wi, % = 1, . . . , N. The structure of these equations resembles the logistic map and they are 
coupled through the average value w(t). The dynamics includes autocatalysis both at the 
individual level and at the community level as well as a saturation term. To model the 
non- stationary conditions we introduce a time dependent parameter into the saturation 
term in each of these equations. We find that the system components spontaneously evolve 
into a power- law distribution of the form of Eq. (JJ), even in the presence of non-stationary 
external conditions. Furthermore, it is shown analytically that the exponent a depends only 
on the ratio of the constant drift component (social security) and the fluctuating component 
(investments). It is found to be insensitive to variations in the saturation term that describes 
the level of economic activity and varies between periods of prosperity and depression. 

The paper is organized as follows. In Sec. II we present the generalized Lotka-Volterra 
model under non-stationary conditions. Analytical results and predictions are presented in 
Sec. Ill and compared with the results of numerical simulations in Sec. IV. A summary is 
presented in Sec. V. 
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II. THE MODEL 



The generalized Lotka-Volterra system [ 131 14, 15, |18, O, |20| describes the evolution in 



discrete time of N dynamic variables tOj, i — 1, . . . , N. In ecological systems, Wi represents 
the population size of the zth specie, while in economic systems it may represent the wealth 
of an individual investor or the market value of a publicly traded firm. At each time step t, 
an integer i is chosen randomly in the range 1 < i < N, which is the index of the dynamic 
variable Wi to be updated at that time step. A random multiplicative factor X(t) is then 
drawn from a given distribution n(A), which is independent of i and t. It will later be 
convenient to express this multiplicative factor by 

X(t) = (A) + n(t) (2) 
where (A) is the average value of 11(A) and 

D = (A 2 ) - (A) 2 (3) 
is its standard deviation. The system is then updated according to 

N N N 

Wi(t +1) = [1 + \(t)]wi(t) + a i,jWj{t) - Y a jti Wi(t) - ^ CijWi{t)wj{t) 

j=i j=i j=i 

Wj (t+1) = Wj (t), j = l,...,N; j^i. (4) 

where a^- and cy are constants. This is an asynchronous update mechanism. The first 
term on the right hand side of Eq. (^), describes the effect of stochastic auto-catalysis at the 
individual level. In an ecological system this term represents variations in the population 
of a given specie, including births and deaths that may be affected by external conditions 
but are not affected by the interaction with other species. In a stock-market system it 
represents the increase (or decrease) by a random factor A(t) of the capital of the investor i 
between time t and t + 1. The second and third terms in Eq. (Q), describe the interaction 
between different dynamical variables. In an ecological system, the second term represents 
the dependence of population i on the availability of food, in the form of population j. The 
third term represents the fact that population i itself may be the food of some other species. 
In an economic system the second and third terms represent trade between investors or 
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firms i and j, such as buying and selling, respectively. The fourth term in Eq. (|j), describes 
saturation effects due to the competition for limited resources. In an ecological model, this 
term implies that large populations tend to exhaust the available resources on which they 
depend. The saturation parameters Cjj are large for populations i and j that consume the 
same type of food. In an economic system this term has to do with the saturation due to 
the finite size of the economy. 

To simplify the analysis we will consider in this paper a simple case in which the W{S 
interact in a uniform fashion with each other. This case is obtained by choosing a it j = a/N 
and Qj = c/N. With this choice Eq. (£|) will be reduced to 



Wi (t + 1) = [1 + X(t)]wi(t) + aw(t) - cwi(t)w(t) 

Wj (t + 1) = wj(t), j = l,...,N; j^i. (5) 

where 

1 N 

=*?£«><(*)■ ( 6 ) 

i=l 

is the average value of the dynamical variables at time t. Here the random term A(i) was 
shifted to \(t) — a but its distribution around the average (and thus the value of the standard 
deviation D) remain unchanged. The second term in Eq. ([5]), may now describe the effect 
of auto-catalysis at the community level. In an economic model, this term can be related 
to the social security policy or to general publicly funded services which every individual 
receives. It prevents an individual u>, from falling below a certain fraction of the average 
w. The third term in Eq. (|5|), describes saturation or the competition for limited resources. 
It has the effect of limiting the growth to values sustainable for the current conditions and 
resources. Within the ecological context, here the interactions between populations are 
uniform, describing the case in which all of them consume the same type of food. We refer 
to Eq. (H]) as the generalized Lotka-Volterra system because when averaged over i and over 
Mt), this system tends to approach a Lotka-Volterra-like equation [p], |2 



w(t + 1) = (1 + (A) + a)w(t) - cw 2 (t). (7) 

where w(t) = w(t). Computer simulations show that after some equilibration time the 
system described by Eq. (H) approaches steady-state conditions. Even at steady state, w 
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exhibits fluctuations. However, its average over long time scales approaches a constant value 
given by 



(w) t = ((\) + a)/c. (8) 

In previous studies of the system described by Eq. ((^) the parameters a and c were 
considered as constants, corresponding to steady conditions of the market. In fact, the 



typical dynamics of microscopic market models [L3|, [TJ], |23j is generically not in a steady 
state. The effect of varying market conditions can be studied by considering the parameters a 
and c and the distribution 11(A) as slowly varying functions of time. We will show below that 
systems described by Eq. ([5]) lead, under very general conditions, to a power-law distribution 
of the Wj's of the form of Eq. (|I]). Moreover, it will be shown that the exponent a is insensitive 
to variations in the parameter c, namely it depends only on a and D. In order to examine 
the effect of variations in the economic conditions we will now introduce an explicit time 
dependence into the third term, as well as a more general dependence on the Wj's. The 
dynamic equation will now take the form 



Wi(t + 1) = [l + \(t))w i (t) + aw(t)-C(w 1 ,...,w n ,t)w i {t) 

Wj(t + 1) = Wj (t), j = l,...,N; j^i, (9) 

where C(wi, . . . , w n , t) is a general function of the w/s that includes an explicit time depen- 
dence. 



III. THEORETICAL ANALYSIS 



In order to study the dynamics of the generalized Lotka-Volterra model, it will be con- 
venient to denote the change of Wj in a single time step by Awi(t) = Wi(t + 1) — w^t). We 
introduce a set of normalized variables 

in ■ 

Xi = —, i = l,..., N. (10) 

w 

The change Axi(t) = Xi(t + 1) — Xi(t) in a single time step is given by 



up to first order in powers of the Aw{8. Considering the time dependence of the average 
w one should remember that at any time step t of Eq. ([|) only one of the w^s is chosen 
randomly and updated. Moreover, there is no correlation between the chosen Wi and A(£). 
Thus, the time evolution of w should be considered on a longer time scale of the order 
of iV moves. However, for simplicity we evaluate Aw by averaging Eq. (|9|) for Awi over 
i = 1, . . . , N at a given time t. We make an independent random choice of X(t) for each i, 
which we denote by Aj(i) = (A) + r)i(t). The time dependence of w is given by 



1 N 

Alu = lyJ2Vj(t)wj(t) + [(A) + a]w(t) - C(w u . . . ,w n ,t)w(t). 



The dynamics of the x^'s is thus given by 



Axi = Xi(t) 



1 



j 



+ a. 



Consider the sum 



TV 



If the Xj f 8 exhibit a distribution of the form 



(12) 



(13) 



(14) 



P(x) ~ x 



-l-a 



then the second moment of the distribution of r(N) satisfies p4|, ^5 



(15) 



(r 2 (N)) 



1/2 



N 1 / 2 , 2 < a 

M 3 -")/ 2 , 1< a < 2 (16) 
N, < a < 1. 

In the first case, the distribution of the x/s exhibits a finite second moment and r(N)/N — > 
in the limit — > oo. In the second case the second moment of P(x) diverges and r(iV) 
follows a Levy distribution. 

In both cases, namely for a > 1, we obtain that in the (thermodynamic) limit — > oo: 



1 

N' 



-x 



(17) 
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Thus, under the assumption that P(x) follows Eq. ([TJJ) with a > 1, we obtain to a good 
approximation that for large values of N 



Axi = [rj(t) - a] Xi + a, i = l,...,N. (18) 

We see that the dynamics of the normalized variable Xi is reduced to a set of identical 
decoupled linear Langevin equations, which do not depend on the function C(wi, . . . ,w n ,t) 
or on the mean value (A) of the multiplicative noise. These equations can be cast into a 
general framework of multiplicative processes of the form 

Ax(t) = r](t)G(x(t)) + F(x(t)). (19) 

Eq. ([TJ) can then be recovered by taking F(xi) = a(l — Xi) and G(xi) = X{. By using a 
suitable change of variables to y — y(x) that satisfies 

^ = — (20) 
dx G{x) 1 ; 

one can reduce Eq. (|19"D to a Langevin equation in which the term rj(t) appears as an additive 
noise, rather than a multiplicative noise such as rj(t)G(x(t)) |26| . The time evolution oiy(t) 
is obtained from Eq. (|T9|) by using the chain differential rule up to second order in Ax (and 
first order in D) 

Ay * d £ Ax + 50 Ai2 - (21) 

Inserting Ax(t) from Eq. (fi9|) and using the change of variables described in Eq. (|2~0D we 
obtain 

Ay ^ ±(F + TfG) + ~(^){F + ^Gf. (22) 

We now approximate the second order term by averaging over the noise term 7], that satisfies 
(rj) = and (rf) = D. We obtain 

F ldG,„ F\ 
A ^" + G-2^< fl+ ^ (23) 

Assuming that F 2 /G 2 D, Eq. (El) is reduced to a discrete-time Langevin equation 
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Ay ~r] + J(y) 
where the drift force J(y) takes the form 



(24) 



Ay) 



F DdG 
G ~ ~2~dx 

The Fokker-Planck equation corresponding to Eq. 



(25) 



is 



dP(y,t) d D d 2 P(y,t) 



(26) 



where P(y,t) is the probability distribution of y at time t. The solution of this equation 
under the stationary condition dP(y,t)/dt = is 



P(y) = exp 



2 rv 



D 



J{y')dyi 



(27) 



Thus, the distribution P{x) = P(y)dy/dx of the original variables is 



P(x) 



G 2 (x 



exp 



2 r Fl - X '\ dx ' 



D J G 2 (x'^ 



(28) 



By taking F(x) = a(l — x) and G(x) = x we obtain a power-law distribution of the form 



P(x) = x 1 a exp 



1 — a 



x 



(29) 



where 



2a 



(30) 



The exponent a thus depends on a single parameter a/D, namely on the ratio between the 
global drift coefficient a and the fluctuations measured by D. 

Another way to derive Eq. (30) from dynamical models of the form (|18|) was shown in 
Refs. H |16 . It is based on the fact that, under steady-state conditions, linear Langevin 
equations of the form (|18D satisfy 0, [16, 29 1 



(fo-o + l) a ) = l. 



(31) 



S 



Considering 77 — a as a small parameter and expanding Eq. ( pi]) in a power series up to 
second order we obtain 



+ < 32 > 



Assuming that a 2 <C D we reproduce Eq. fl30|). 



IV. NUMERICAL SIMULATIONS AND RESULTS 

To examine the theoretical predictions presented in Sec. Ill we have performed computer 
simulations of the generalized Lotka-Volterra system described by Eq. (|9|) with different 
choices of C(w\, . . . , w n , t). It was found that after some equilibration time the distribution 
P(x) reaches a steady state, and exhibits a power-law behavior. For C(w\, . . . ,w n ,t) = cw 
[as in Eq. @], w(t) fluctuates around some average value, given by Eq. (||). For a general 
function C(w\, . . . , w n , t), that exhibits an explicit time dependence, w(t) continues to vary 
according to this function and its temporal average does not reach a steady state. 

To examine how robust the power-law distribution is under varying conditions we have 
simulated Eq. @ with 

/ 27rt\ N 

C(wi, . ..,w n ,t) = c \1 + sin ^rj ^j^J ( 33 ) 

(for Co = 0.001, T = 2 x 10 5 and (A) = 0.002) and compared the results with the case 
C(wx, . . . ,w n ,t) = cw (for c = 1 and (A) = 0.01). The time dependence of w in both 
cases is shown in Fig. 0(a). The distributions P(x) obtained from the simulations in these 
two cases are shown in Fig. 0(b). The distributions are found to be nearly identical and 
exhibit a power-law behavior characterized by the same exponent a. The exponent a is also 
found to be independent of (A). Note that the power- law behavior is maintained even for 
C(wi, . . . ,w n ,t) = 0, where w(t) does not reach a steady state and diverges to infinity (for 
(A) > 0) or collapses to (for (A) < 0) [[TJJ. This can lead to changes by orders of magnitude 
in the total wealth or the population size without affecting the exponent a. 

To examine the theoretical prediction for the distribution, given by Eq. (^), and the ex- 
ponent a, given by Eq. ([30]) , we have compared these predictions to the results of numerical 
simulations. This comparison for the distribution P(x) is shown in Fig. |2|. The simulations 
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were done for N = 1000, a = 0.00023, C(wi, . . . , w n , t) = O.Olw), and A(t) uniformly dis- 
tributed in the range 0.0 < A(t) < 0.1 (namely, D = 0.00083). We found that the simulation 
results are in very good agreement with the theoretical predictions. A power law distribution 
is found for a range of almost three orders of magnitude, with the exponent a = 1.52. This 
is close to the theoretical prediction of a = 1 + 2a/ D = 1.55. The distribution has a peak 
at Xq = (a — l)/(a + 1), that using Eq. ( |3"0"D can be expressed by xq = a/(a + D). Above 
Xo the distribution P{x) behaves like power-law while below it P(x) decays exponentially. 
This provides an effective lower cutoff for the range of x in which a power-law behavior 
is observed. This result can be compared to a somewhat simpler model studied earlier, in 



which the value of the lower cutoff x m i n is imposed as a constraint [15|. In this model, using 
the sum rules for the probability and the total wealth, it was found that x m i n = 1 — 1/a. 
Using Eq. (|30"D it can be expressed as x m i n = 2a /{2a + D). These predictions for the lower 
bounds in the two models satisfy Xq < x m i n < 2xq, namely, they are in good agreement in 
light of the broad distribution of x. 

To examine the prediction given by Eq. (|30| ) for the exponent a, we present in Fig. |3| a 
comparison between this prediction and the numerical results for a as a function of a/D. 
The numerical results are presented for iV = 1000 (•). The prediction of Eq. ( pOf ) (solid line), 
shows a good agreement with the numerical results for a/D > 0.2. The numerical results for 
the range of small a/D converge to the theoretical prediction as the value of N is increased. 
As shown in Ref. \TE\ the infinite system limit, N —>■ oo, and the vanishing coupling limit, 
a/D — > 0, do not commute. On the one hand, for any finite iV and a/D — ► the exponent 
a — > 0. On the other hand, for any fixed positive value of a/D (no matter how small) and 
iV — >• oo the exponent a > 1. The majority of empirical results in which 1 < a < 2, indicate 
that the second case is highly relevant and that the theoretical predictions of Eqs. (|29| ) and 
(j30|) broadly apply. Thus, for iV — > oo, both the exponent a of the power-law decay, and 
the lower bound xq depend only on a single parameter a/D. In the economic context, this 
parameter represents the ratio between the fixed income of minimal wage jobs or social 
security payments and the level of fluctuations of the speculative income/loss. 
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V. SUMMARY 



We have studied the dynamics of stochastic Lotka-Volterra systems under non-stationary 
conditions using both analytical and numerical techniques. For this class of models, we found 
that in order to obtain a power-law distribution, it is sufficient that relative returns of the 
agents are stochastically equivalent. The assumption that the distribution 11(A) of the mul- 
tiplicative noise, is independent of i, means that there are no investors or strategies that can 
obtain 'abnormal' returns. This can be related to to the 'efficient market hypothesis', which 
assumes that the market pricing mechanism is so efficient that it reaches the 'right price' be- 
fore any of the agents can take systematic advantage. Therefore, the presence of a power-law 
distribution may be a sign of 'market efficiency', by analogy with Boltzmann distributions in 
statistical mechanics systems, which characterize thermal equilibrium. Here we have shown 
that the power-law distribution is stable even under non-stationary economic conditions, 
that are represented by the time dependence of the saturation term C(w\, . . . ,w n ,t). We 
found that even under such conditions the distribution of the (normalized) dynamical vari- 
ables Xi follow a power-law distribution with an exponent a. An expression for a in terms 
of ratio of the parameters a and D was obtained [Eq. ([30])]. In the economic context, the 
parameter a represents the minimal wage or social security payments, while D represents 
the level of fluctuations in speculative income/loss. These results provide the distribution 
of wealth in a society in terms of the social security policy and the volatility of the stock 
market. They also provide a connection between the incomes/wealths of the poorest and 
the richest sectors of the society as a function of a single parameter. 



11 



[1] V. Pareto, Cours d'Economique Politique, (Macmillan, Paris, 1897), Vol 2. 

[2] B. Mandelbrot, Econometrica 29, 517 (1961). 

[3] B. B. Mandelbrot, Comptes Randus 232, 1638 (1951). 

[4] B. B. Mandelbrot, J. Business 36, 394 (1963). 

[5] A.B. Atkinson and A.J. Harrison, Distribution of Total Wealth in Britain (Cambridge Uni- 
versity Press, Cambridge, 1978). 
[6] H. Takayasu, A.-H. Sato and M. Takayasu, Phys. Rev. Lett. 79, 966 (1997). 
[7] P. Jogi, D. Sornette and M. Blank, Phys. Rev. E 57, 120 (1998). 
[8] M. Levy and S. Solomon, Physica A 242, 90 (1997). 
[9] R.N. Mantegna and H. E. Stanley, Nature 376, 46 (1995). 
[10] P. Gopikrishnan, M. Mayer, L.A.N Amaral and H.E. Stanley, Euro. Phys. J. B 3, 139 (1998). 
[11] V. Plerou, P. Gopikrishnan, B. Rosenow, L.A.N Amaral and H.E. Stanley, Phys. Rev. Lett. 
83, 1471 (1999). 

[12] P. Gopikrishnan, V. Plerou., L.A.N. Amaral, M. Meyer and H.E. Stanley, Phys. Rev. E 60, 
5305 (1999). 

[13] M. Levy, H. Levy and S. Solomon, J. Phys. 5, 1087 (1995). 

[14] S. Solomon and M. Levy, Int. J. Mod. Phys. C 7, 745 (1996). 

[15] O. Malcai, O. Biham and S. Solomon, Phys. Rev. E 60, 1299 (1999). 

[16] D. Sornette and A. Johansen, Physica A 245, 411 (1997). 

[17] D.H. Zanette and S.C. Manrubia, Phys. Rev. Lett. 79, 523 (1997). 

[18] O. Biham, O. Malcai, M. Levy and S. Solomon, Phys. Rev. E 58, 1352 (1998). 

[19] M. Levy and S. Solomon, Int. J. Mod. Phys. C 7, 595 (1996). 

[20] O. Biham, Z.-F. Huang, O. Malcai and S. Solomon, Phys. Rev. E 64, 026101 (2001). 

[21] Elements of physical Biology, edited by A.J. Lotka (Williams and Wilkins, Baltimore, 1925). 

[22] V. Volterra, Nature 118, 558 (1926). 

[23] M. Levy, H. Levy and S. Solomon, Europhys. Lett. 45, 103 (1994). 
[24] M.F. Shlesinger, G.M. Zaslavsky and J. Klafter, Nature 363, 31 (1993). 
[25] J. Klafter, G. Zumofen and M.F. Shlesinger, in Levy flights and related topics in physics, 
edited by M.F. Shlesinger, G.M. Zaslavsky and U. Frisch (Springer- Verlag, Berlin, 1995), pp. 



12 



196-215. 

[26] P. Richmond, Euro. Phys. J. B 20, 523 (2001). 

[27] N.G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amster- 
dam, 1981). 

[28] S. Solomon and P. Richmond, Physica A 299, 188 (2001). 
[29] H. Kesten, Acta. Math. 131, 207 (1973). 



13 



FIG. 1: (a) The time dependence of the average wealth w for the model of Eq. @ with 
C{w\, . . . , w n , t) given by Eq. fl33|) where Co = 0.001 and T = 2 x 10 5 (upper curve), and with 
C(wi, . . . ,w n ,t) = cw where c = 1 (lower curve). In the first case w oscillates, following the time 
dependence of C(w\, . . . , w n , t), while in the second case it only exhibits small fluctuations around 
a constant value. In both cases N = 100, a = 0.00083, D = 0.0033; (b) The distributions of 
the variables Xi = Wi/w for the simulations shown in the lower curve (squares) and the upper 
curve (•) in (a). The two distributions are found to be nearly identical, showing an approximate 
power-law behavior. We thus observe that the exponent a is robust and insensitive to variations 
in C(wx, ■ ■ -,w n ,t). 

FIG. 2: Results of computer simulations (dots) and theoretical analysis based on Eq. (^9) (•) 
for the distribution of the variables Xi = Wi/w. The parameters are N = 1000, a = 0.00023, 
D = 0.00083 and C(wi, . . . , w n , t) = O.Olu). In both cases a power-law distribution is obtained 
with an excellent agreement between the theoretical predictions and the simulation results. 

FIG. 3: Simulation results for the exponent a of the power-law distribution of the variables 
Xi = Wi/w, i = 1, . . . , N as a function of the parameter a/D for N = 1000 and C(w\, . . . , w n , t) = 
O.OOOOlw (•). The theoretical prediction of Eq. (^) (line) is found to be in agreement with the 
numerical values for a/D > 0.2. The agreement for small values of a/D tends to improve as N 
increases. 
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